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(57) Abstract: The wavefront shape 
of radiation arriving at a pupil plane 
(4) is defined by a set of predetermined 
orthonormal functions such a Zernike 
functions with a weighting coefficients 
as a function of the shape. Apparatus 
for determining at least one such 
weighting coefficient comprises rate 
means responsive to radiation as 
it traverses the input pupil (4) to 
determine a pixelwise distribution 
indicative of rate of radiation intensity 
change thereat, and converting means 
for converting the intensity distribution 
to one or more weighting coefficients. 
The converting means comprises a 
store holding one or more matrices 
of predetermined values, each said 
matrix corresponding to one said 
orthonormal function, and the size 
of each said matrix corresponding to 
the number of pixels in the pixelwise 
distribution, and calculating means for 
multiplying said pixelwise distribution 
by a said matrix and adding the results 
to provide the weighting coefficient for 
its orthonormal function. The rate of 

radiation intensity change at pupil (4) is conveniently approximated by measurements of intensity in closely adjacent planes (22, 
23) either side of pupil (4). 
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Wavefront Sensing 

The present invention relates to the determination of properties of a wavefront as it 
reaches an input pupil, for example of an imaging system or a measuring system. 

Our copending International Patent Application No. WO 99/46768 (published in the 
5 name of the Secretary of State for Defence, inventors Paul Blanchard and Alan 
Greenaway) describes a 3 -dimensional imaging system utilising a diffraction grating 
which is distorted according to a quadratic function. When located adjacent an 
imaging system such as a biconvex lens, the zero order light from the grating provides 
the normal focussed image in an imaging plane via the imaging system, and the 
10 positive and negative diffraction orders from the grating, particularly the first orders 
(although higher orders could contribute) provide respective focussed images in the 
same imaging plane but laterally displaced by respective amounts in opposed 
directions from the normal image. 

A related article by the same inventors, "Simultaneous Multi-plane Imaging with a 
15 Distorted Diffraction Grating" Applied Optics 38 (32), 6692-6699 (1999) also 
mentions applications of such a system, including the improvement of imaging 
systems by sensing phase aberration introduced for example by the atmosphere and 
using corrective adaptive optics to improve the focussed image. It will be understood 
that atmospheric and other effects will lead to a distortion of a wavefront as it 
20 propagates, and that this includes a variation in phase with solid angle over that 
ideally expected. 

More generally, it is desirable to be able to determine wavefront shape (which may be 
alternatively regarded as distribution of local phase), either as an absolute 
measurement or as an aberration from an ideal shape or phase distribution. Many 

25 known methods are based ofrthe a posteriori reconstruction of a complex amplitude, 
from energy measurements or from phase-unstable data. Such reconstructions 
commonly are complex and time consuming. However, the need in some fields for 
fast accurate and robust determinations of wavefront shape or phase aberrations has 
led to a number of real time techniques. This is so particularly but not exclusively in 

3 0 the optical field, with the increasing application of adaptive optics (AO), using real- 
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time wavefront modulation to correct stochastic aberrations (usually, but not 
necessarily, turbulence induced). 

At least for some applications such as the aforementioned adaptive optical system, it 
is even more desirable to be able to do so in what is essentially real time. One such 
5 technique involves the Shack-Hartmann sensor in which the wavefront is passed 
through an array of lenslets each of which forms a sub-image of the scene. The 
displacement of the sub-images from the local axial position for each lenslet in the 
array provides a matrix of local, two-dimensional wavefront gradients, from which 
the wavefront can be reconstructed by integration. This method has been widely and 
1 0 successfully used in astronomy. 

However, if the scene whose image is to be corrected for aberrations is extended and 
low-contrast, the detector format for each sub-image, and the computational effort 
required to compute the displacements of the sub-images from their local axes 
(usually by cross correlation of the sub-images), can become prohibitive. Other 
15 approaches, such as image sharpness and wavefront-shearing interferometry have 
been used to circumvent these problems with varying degrees of success. 

It is known that the shape of a wavefront can be described by a set of appropriately 
weighted orthonormal functions. The functions themselves are predetermined, but the 
weighting is variable, and it is the weighting coefficients which define the shape. 

20 It is also known that the wavefront shape may be determined by application of the 
intensity transport equation. This can be a complex and difficult process, which is 
dealt with in embodiments of the present invention by using certain arrangements for 
detecting the radiation, and by providing for each weighting coefficient a 
predetermined matrix of values for use in real time measurement of the coefficient. 

25 The matrix of values is related or matched to the arrangement for detecting the 
radiation, and determination of the matrix, which remains a complex and lengthy 
process, is effected initially, i.e. prior to the radiation detection. 

The present invention relates to a version of an established algorithm known as 
wavefront curvature sensing, which may be regarded as a special case of phase 
3 0 diversity (see, for example, Gonsalves R A, "Phase retrieval and diversity in adaptive 
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optics", Opt. Eng. 21 (1982) 829-832; and Roddier F, Curvature sensing and 
compensation: a new concept in adaptive optics, Appl. Opt, 27 (1988) 1223-1225, 

Generally, 'phase-diverse' solutions to the problem of wavefront reconstruction from 
two intensity measurements are related to the 'two-defocus* algorithm described in 
5 Misell D L, "An examination of an iterative method for the solution of the phase 
problem in optics and electron optics", J. Phys. D.; Appl. Phys. 6 (1973) 2200-2216. 

In a first aspect the present invention provides a measurement method for determining 
data relating to the local shape (or distribution of local phase) of a radiation wavefront 
arriving at a pupil plane, wherein said shape is defined by a set of predetermined 

10 orthonormal functions, each function being provided with a weighting coefficient for 
determining the shape, said data comprising at least one said weighting coefficient, 
the method comprising determining a pixelwise distribution indicative of rate of 
change of radiation intensity as the radiation traverses the input pupil, and converting 
said pixelwise distribution to said data, wherein said converting step comprises 

15 providing one or more matrices of predetermined values, each said matrix 
corresponding to one said orthonormal function, and the size of each said matrix 
corresponding to the number of pixels in said pixelwise distribution, for each said 
matrix multiplying said pixelwise distribution by said matrix and adding the result to 
provide said weighting coefficient for its said orthonormal function. 

20 The first aspect of the invention extends to measuring apparatus for determining data 
relating to the local shape (or distribution of local phase) of a radiation wavefront 
arriving at a pupil plane, wherein said shape is defined by a set of predetermined 
orthonormal functions, each function being provided with a weighting coefficient for 
determining the shape, said data comprising at least one said weighting coefficient, 

25 the apparatus comprising a said input pupil, means for determining a pixelwise 
distribution indicative of rate of change of radiation intensity as the radiation traverses 
the input pupil, and converting means for converting said intensity distribution to said 
data, wherein said converting means comprises a store holding one or more matrices 
of predetermined values, each said matrix corresponding to one said orthonormal 

30 function,, and the size of each said matrix corresponding to the number of pixels in 
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each of said first and second images, and calculating means for multiplying said 
pixelwise distribution by a said matrix and adding the results to provide said 
weighting coefficient for its said orthonormal function. 

As will be described in a largely theoretical section later, it is necessary to know how 
5 the local intensity of the radiation across the input pupil is changing as it passes 
through the input pupil. As particularly described, this is not effected directly, but by 
measuring the distributions of intensity in respective secondary planes either side of, 
and adjacent to, the input pupil, and measuring the distribution of intensity differences 
between the two distributions. 

10 The intensity distributions in the two secondary planes may be may be obtained by 
any method known or obvious to the skilled person, for example by the use of a 
conventional beam splitter placed after the input pupil and directing light to respective 
detectors or detector arrays. This is shown schematically in Figure 5 in relating to 
spaced planes, a beam splitter 24, a ccd camera 26 for receiving a focussed image of 

15 plane 22 via a lens 25, and a ccd camera 28 for receiving a focussed image of plane 
23 via a lens 27. A possible difficulty is synchronisation of the two cameras, and of 
course there is the need to* ensure and to maintain alignment and calibration. 

However, in a particularly preferred optical embodiment the two intensity 
distributions are provided on a common plane. For example in an optical system a 

2 0 diffraction grating of the general type described in our copending International Patent 
Application WO 99/46768 mentioned above may be used in conjunction with an 
imaging lens or other optical imaging system for directing and focussing the light 
from at least the plus one and minus one orders, and optionally the zero order, onto 
the plane of a single detector or detector array. The focussed light received by the 

25 detector (array) includes simultaneous images of the laterally displaced images of the 
two secondary planes (plus and minus one orders), and optionally the zero order. The 
former images are employed to provide the aforesaid distribution of intensity 
differences. 
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A computer generated hologram may be employed instead of the aforesaid grating for 
the same purpose, and with or without an additional lens according to the properties 
of the hologram (which can itself additionally provide the lens function). 

Other arrangements for obtaining both intensity distributions in a common plane are 
5 described in our copending UK Patent Application No. GB 0301923.9 (ref: 
P21386GB). 

Thus in a preferred embodiment the said pixelwise distribution is obtained by (means 
for) providing focussed first and second images of respective first and second 
secondary planes lying adjacent, and respectively before and after, the pupil plane, 
10 and deriving the pixelwise distribution as the pixelwise distribution of intensity 
difference between the first and second images. The first and second images may 
alternatively be regarded as images of the pupil plane with relatively minor degrees of 
positive and negative defocus. 

The wave is not necessarily an optical wave, and other exemplary areas where a 
1 5 knowledge of phase aberrations may be useful include studies on x-ray diffraction and 
nuclear structure (see for example von Kampen N G, "S-matrix and causality 
condition I. Maxwell field", Phys. Rev. 89(1953)1072-9); subatomic particles such as 
in electron microscopy (see for example Misell D L, "An examination of an iterative 
method for the solution of the phase problem in optics and electron optics", J. Phys. 
20 D: Appl. Phys. 6(1973)2200-2216; and Gerchberg R W and Saxton W O, "A practical 
algorithm for the determination of phase from image and diffraction plane pictures", 
Optik 35(1972) 237-246); and phase-unstable interferometry (see for example 
Schwarz U.J., "Mathematical-statistical description of the iterative beam removal 
technique (method CLEAN)", Astron. Astrophys. 65(1978)345-356; and Lannes A, 

2 5 "Backprojection mechanisms in phase-closure imaging", Exp. Astron. 1(1989)47-76). 

As particularly described and as implemented at present the radiation is optical (ultra- 
violet to infra-red), and preferably visible or infra-red. 

The set of orthononnal functions may be Zernike modes, although other sets could be 
employed, depending inter alia on the potential application of the data so obtained and 

3 0 the nature of the radiation. 
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For later processing, it is necessary to provide the distribution of intensity differences 
in the form of a spatially discrete array of values. Accordingly, it is preferred to 
employ a pixellated detector array (although "pixel" in the present context refers to an 
image of the input pupil as opposed to the more normal reference to an object plane). 
5 However, as will be recognised by the reader, it is possible to use a non-pixellated 
detector, e.g. a conventional CRT camera with subsequent image processing to 
provide a like result. References herein to a detector should be taken to embrace both 
spatially pixellated (or digitised) detectors, as well as those providing a continuous 
indication of intensity with position. 

10 The aforementioned use of beam splitters is not only costly, but can give rise to 
problems in realising the necessary precision in alignment, particularly insofar as it is 
necessary to ensure accurate relative alignment of the two near-pupil plane intensity 
distributions. The use of the distorted diffraction grating or the hologram can not only 
lead to lower costs, but it can help avoid any relative misalignment between the two 

1 5 images of the secondary planes. 

The number of weighting coefficients provided in use of the method according to the 
invention can vary from one upwards. If the discrete or "pixellated" array of intensity 
difference comprises N values, a corresponding number N of weighting coefficients 
for N respective orthonormal functions can be obtained if required. While this may 

20 involve a relatively high number of individual computational steps, they are simple 
steps capable of being effected quickly. For example, with a 40 by 40 array, giving 
1600 intensity difference values effectively across the input pupil, there will be 1600 
matrices for 1600 orthonormal functions, each function requiring 1600 
multiplications of its matrix values by the 1600 intensity difference values, followed 

25 by an addition step. 

It should be clear that the fineness of the structure of the phase distribution to which 
the method of the invention relates is at least in part determined by the number of 
values in the discrete or "pixellated" array of intensity difference. However, in many 
applications it is not necessary to determine all of the possible orthonormal weighting 
3 0 coefficients. For example, in the above mentioned use in correcting optical images 
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for atmospheric distortion a corrective or adaptive light modulator placed in the 
optical imaging path may only be capable of dealing with a limited number of modes 
corresponding to a like number of the orthonormal function, so that it is only 
necessary to determine the corresponding limited number of weighting coefficients. 
5 In one particular embodiment, the light modulator comprises a distortable membrane 
mirror capable of phase correcting for 37 modes, and so it is only necessary to 
determine the corresponding 37 weighting coefficients. If another light modulator is 
employed for phase correction, for example a micro-electro-mechanical (MEM 
displaceable mirror array, e.g. of 8 by 8 elements, the number of correctable modes 
1 0 might correspond to the number of elements, i.e. 64. 

Only one coefficient may be required in certain cases, such as in range finding or the 
determination of the degree of focus (see later). 

The first two Zernike modes are tip and tilt. To correct for these modes, e.g. to 
correct the pointing of an imaging system or a laser beam, a simple tip/tilt movable 
15 mirror can be employed. Other Zernike modes will require more complex 
arrangements such as those just exemplified. 

Another reason for not determining the weighting coefficients for a larger set of the 
orthonormal functions may be the time periods involved. As the phase structure 
becomes more finely resolved, the fine structure progressively relates to spatially 
20 small features, whereas a lower number of functions would relate to larger scale 
features. It is to be expected, or likely, that large scale features change slowly relative 
to the smaller scale features. In a real time application such as image correction, there 
is little point in determining the weighting coefficient for a function if that weight 
coefficient is expected to change at a rate comparable to or faster than the rate at 

2 5 which the coefficient can be determined. 

While the measurement of the weighting coefficient(s) may be useful in itself, for 
example in giving a measure of the degree of phase perturbation, or the magnitude of 
the cause giving rise to such perturbation, one particular application of the method is 
to imaging, and in particular to correcting for phase perturbations arising as the 

3 0 radiation travels towards the input pupil of an imaging system. For example, visible 
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light is commonly affected by atmospheric turbulence and density changes, 
particularly when the viewed object is relatively distant, and further perturbation may 
arise for example if scintillation occurs, or there is atmospheric pollution. As 
indicated previously, the focussed image can be improved by determining phase 
5 aberrations in the received wavefront and using corrective adaptive optics in the 
imaging system, but this is difficult in real time. Thus one reason for the development 
of the present invention is to meet a need for real-time wavefront reconstruction in 
terrestrial imaging applications, where the scene under observation is extended and 
intrinsically low contrast and where the use of laser beacons to give a wavefront 
1 0 sensing reference is prohibited because of eye safety or other issues. 

The corrective optics means may be used independently of the phase measurement 
(for example by using a beam splitter), but in a preferred embodiment the wavefront 
measuring system and the imaging system receive light from the same the corrective 
optics means, and the arrangement is such as to minimise the measured phase 
15 aberrations. 

Thus the invention extends to an imaging method (preferably, but not necessarily, 
optical, and most preferably visible or infra-red imaging) in which weighting 
coefficient(s) are obtained according to the first aspect of the invention, and wherein 
corrective optics are incorporated in the imaging system and are controlled in 

20 response to the determined weighting coefficients. The corrective optics can be a 
spatial phase modulator, for example a deformable mirror membrane or a micro- 
mechanical movable mirror array, optionally with a tip/tilt mirror for dealing 
separately with the first two Zernike modes. These modes have the potential for 
involving relatively large scale angular corrections, and a separate mirror, is therefore 

25 preferred, although, if necessary, this mode could be dealt with by the same spatial 
phase modulator as deals with other mode(s). 

In the simplest implementation, the corrective optics may be peculiar to the path for 
the imaging radiation. That is to say, the weighting coefficients are determined using 
light which does not encounter the phase control means. 
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However, in a preferred embodiment, the phase control means is common to the 
imaging radiation and that used for determination of the weighting coefficients. 
Operation of this embodiment maybe such that the weighting coefficients are reduced 
as far as possible, and preferably completely nulled, by operation of the phase control 
5 means in a feed back loop. It will be appreciated that the signal in the feed back loop 
remains indicative of the weighting coefficients determined according to the first 
aspect of the invention. 

In a similar way to the above application of adaptive imaging, a light beam, for 
example a laser beam, may be controlled. The beam wavefront may be measured 
10 according to the invention and adaptive optics used to control it. As in the imaging 
system the adaptive optics lies in a main optical path, and the measuring system may 
be located either in the main path e.g. in a null sensing system, or in a separate path 
(for example from a beam splitter). 

Another application of the invention is to the determination of atmospheric 
15 turbulence, which is a function of the measured coefficients, and which may be used 
for example in the correction or control of other (non-imaging) optical measurements. 

A further application is in the field of focus or range determination and autofocus. 
When the wavefront mode being sensed is the component of defocus a (Zernike 
polynomial Z 2 ,o(r,9)), the range of the source is related to the coefficient of defocus a, 
20 by 



R' = R 



ttD 2 



(0.1) 



where R is the front focal length of the system (i.e. the range at which an object is 
perfectly focussed), and D is the diameter of the pupil. 

The invention works well when the viewed field contains a single high intensity 
25 point, enabling well defined and spatially separate images of light from the point as in 
the secondary planes either side of the pupil plane to be formed on a single sensor 
surface. If no such high intensity point exists, it is possible to provide one, for 
example by suitable illumination for example using a laser beam. 
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Where there are no such singular points there may be overlap of the images on the 
sensor plane, and the measurement process may be correspondingly difficult. 

In one embodiment, a plurality of points in the viewed field are illuminated, for 
example using a light source or laser with an appropriate diffraction grating. It is 
5 arranged that the corresponding two sets of points on the sensor plane (i.e. from the 
planes either side of the input pupil) are spatially separate. In this way measurements, 
for example of range, can be taken at a number of points in the field of view. 

In the invention, the provision of one or more predetermined matrices facilitates 
essentially real time calculation of the weighting coefficient(s) since the processing 

10 amounts to image intensity subtraction or comparison, matrix multiplication of the 
resulting intensity distribution and addition of the results, all of which are simple 
steps capable of being performed very rapidly on any modern computer. By contrast, 
the calculation of the matrices themselves is a much more complicated process taking 
several hours or longer to perform, and is performed prior to sensing of a radiation 

15 wavefront. The method of the invention has been implemented experimentally and 
has been found to provide fast, accurate and robust wavefront reconstruction. 

Further features and advantages of the invention will become clear upon a perusal of 
the appended claims, to which the reader is referred, and upon a reading of the 
following more detailed description of an embodiment of an imaging system 
2 0 according to the invention including apparatus according to the invention, made with 
reference to the accompanying drawings, in which: 

Figure 1 schematically illustrates an adaptive optics imaging system employing 
apparatus according to the invention; 

Figure 2 schematically illustrates the principle behind the determination of the rate of 
2 5 change of radiation (light) intensity at an input pupil; 

Figure 3 schematically illustrates an arrangement for determining the pixelwise rate of 
change of intensity of radiation as it passes between a pair of closely spaced planes 
either side of an input pupil; 
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Figure 4 schematically illustrates an alternative arrangement for determining the 
pixelwise rate of change of intensity of radiation as it passes between a pair of closely 
spaced planes either side of an input pupil; 

Figure 5 shows an alternative method of obtaining the requisite first and second 
5 focussed images on separate CCD cameras; and 

Figure 6 is an image for illustrating a method of correcting for background radiation. 

In the system illustrated in Figure 1, it is desired to focus an optical image of a distant 
object 10 in an object plane 1 in a focal plane 2, but without correction the focussed 
image tends to be distorted by atmospheric turbulence 3 acting on the light 31 before 
10 it reaches the input pupil 4. 

Light from the input pupil 4 is reflected by a tip/tilt mirror 5 (i.e. a mirror which can 
be moved angularly about two axes), and reflected by a spatial phase modulator 6 
comprising a deformable mirror. Light from modulator 6 is passed through a 
beamsplitter 7 for reflection by a fixed mirror 8 and transmitted by an imaging lens 

15 acting to focus light from the object plane 1 on the image plane 2. Light reflected by 
the splitter 7 is passed to and detected by a wavefront sensor 1 1 which determines the 
weighting coefficients for the first (tip/tilt) and other selected Zernike modes, and 
produces respective signals 12, 13 for correspondingly controlling the tip/tilt mirror 5 
and the modulator 6. The mirror 5 and modulator 6 are common to the optical paths 

20 to the focal plane 2 and the sensor 11, and the latter is arranged to minimise or null 
the weighting coefficients in the wavefront of the sensed light, so that the optical 
wavefront 14 transmitted towards the focal plane is rendered more ideal - that is 
substantially planar, having had phase perturbations or distortion relative to a planar 
wavefront substantially removed by the action of the mirror 5 and modulator 6. 

25 In a variation of the embodiment of Figure 1, the light 31 is a laser beam with 
potential imperfections in the wavefront, not necessarily due to atmospheric 
turbulence. The, or some, imperfections are removed in the output beam as indicated 
at 14. Reflector 8 and lens 9 may not be required. 
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The embodiment of Figure 1 or its variant may be further varied by relocating the 
beam splitter 7 and sensor 1 1 before the tip/tilt mirror 5 and modulator 6 and using 
the output of sensor 11 for forward control of mirror 5 and modulator 6 without the 
use of a nulling system. 

5 In an imaging system of the type shown in Figure 1 or its variant, the Zernike 
polynomial Z2,o(r,9) pertaining to defocus may be employed either to control the 
modulator 6, or to control the position of the lens 9, for auto-focus. This coefficient 
may also or alternatively be employed to provide the range of an object in the field of 
view. 

1 0 Figure 2 schematically illustrates the principle behind the determination of the rate of 
change of radiation (light) intensity at the input pupil 4. This is discussed in more 
detail in the theoretical part of the description to follow, but for now it is noted that 
phase distortion effectively are manifested as distortions from pianarity (in this case, 
but any standard shape, e.g. spherical with a predetermined rate of curvature, could be 

15 used) and that such distortions cause the radiation to converge or diverge, giving rise 
to a variation in local intensity as the radiation travels between two planes 22, 23 
closely equally spaced either side of the pupil plane (not shown). This variation can 
be measured and presented in a pixelwise manner for use in the apparatus of the 
present invention, e.g. in the sensor 11 of Figure 1. 

20 Figure 3 shows how the optical input to the sensor 1 1 could be derived. A focussing 
lens 15, such as the biconvex lens illustrated, is arranged to focus an image of an 
object plane 16 on a measuring plane 17 of an image sensor such as a CCD or other 
digital camera array. Closely adjacent the lens (for example before it, as shown) is 
provided a quadratically distorted diffraction grating 18 of the type generally 

25 described in our copending International Patent Application WO 99/46768. Together 
with the zero order component, which provides an image of the plane 16 centrally on 
plane 17 at a position 19, the grating provides usefixl plus one and minus one 
diffraction orders at least. These orders respectively effectively slightly decrease and 
increase the power of the lens 15, so as to produce in the plane 17 respective focussed 

3 0 images of planes 22 and 23 axially either side of the image plane 16, at locations 20 
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and 21 laterally either side of the position 19. Planes 22 and 23 are located close to 
the plane 16. 

It is thus possible to measure the pixelwise intensity distribution of the light at planes 
22 and 23, and pixelwise to subtract the resulting measurement to obtain an indication 
5 of the pixelwise rate of change of light intensity at the pupil plane 1 6. 

There are disadvantages attached to the use of the arrangement of Figure 3, one 
principal disadvantage being that the magnification of the images at locations 19 to 21 
differ. An alternative and preferred arrangement for providing an optical input to the 
sensor 11 is shown in Figure 4, where the lens 15 is set to focus sources at infinity on 

10 a central focussed spot 19' the sensor 17. The grating 18 acts with the lens 15 so that 
an image of a closer real image plane A is provided by the plus one diffraction order 
at a position 20' on the sensor 17 on one side of the central position 19', and an equal 
magnification image 2V of a virtual plane is provided by the minus one diffraction 
order at an equal distance on the other side of the central position 19'. The virtual 

15 plane is shown as a plane C, shown on the other side of the lens 15, but the skilled 
person will realise that this represents an image plane effectively beyond infinity, i.e. 
the images 20' and 21' effectively relate to planes either side of the pupil plane. The 
images 20' and 21 9 are used in lieu of the images 20 and 21 of Figure 3. 

To provide a somewhat more detailed explanation of the reasoning behind the 
2 0 invention, it is necessary first to consider the Intensity Transport Equation (I.T.E.) 

For light of wavelength X propagating along the z-axis in the paraxial approximation, 
the complex amplitude w 2 (r) satisfies the parabolic equation 



25 The complex amplitude may be expressed in terms of the real-valued intensity, 7 z (r), 
and phase, <|> z (r), as follows: 




(1) 



where V 2 = 
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«,(r) = [/,(r)f exp[/<4,(r)] 

Multiplying Equation (1) on the left by «/(r) and the complex conjugate of Equation 
(1). on the left by w 2 (r), and subtracting the resulting expressions, leads to the 
Intensity Transport Equation (ITE), 

5 -^/,(r) = V.(/ 1 (r)V^(r)) (2) 

oz 

To help in understanding the meaning of the I.T.E., Equation (2) can be written 

-^^/ z (r) = /,(r)V 2 ^(r) + V/ J (r).V^(r) (3) 

oz 

This contains a "curvature" term / z (r)VV z (r), and a "slope" term W z (r)-V^(r). 
In a region where J z (r) is uniform the slope term vanishes (VJ,(r) = 0) and Equation 
10 (3) simplifies to: 

(4) 

I dz 

which relates the local increase or decrease in intensity as light propagates along the 
optical axis to the local curvature of the wavefront, as indicated in Figure 2. Thus, if 
the intensity distribution and its derivative in the z-direction are measured in the pupil 
15 plane of an optical system, the wavefront curvature can be obtained from Equation 
(4). 

In practice, measurement of the intensity in two closely-spaced planes symmetrically 
placed with respect to the entrance pupil provides a practical method for the 
estimation of the axial intensity gradient. The schematic in Figure 2 illustrates this 
20 principle, which is similar to reconstructing the shape of waves on the surface of a 
swimming pool from the scintillation patterns observed on the pool floor. 

Knowledge of the wavefront curvature alone is sufficient to determine the wavefront 
only to within an arbitrary harmonic function. Boundary conditions are also required 
in order to obtain a unique solution. If the region in which we wish to reconstruct the 
2 5 wavefront consists of a uniformly illuminated aperture, we have: 
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where 

I 0 is the constant intensity, 
A denotes the system aperture, 
5 C is the aperture perimeter, 

W A is the aperture function (=1 inside^, =0 outside A\ 
8c is a delta-function around C, 
n is an outward-pointing unit vector normal to C. 
Equation (3) then becomes 

i o ~ k iw = W-Vcfi • V<* z (r) . (5) 

Measurement of the axial intensity derivative thus consists of both wavefront 
curvature within the aperture and wavefront slope around the aperture edge. This 
information is sufficient to determine the phase uniquely to within an arbitrary 
additive constant, which is consistent with an inability to measure absolute phase. 

15 The above equation can be solved using Green's 2 nd identity. Lety(r) and g(r) be any 
two continuous functions of r having finite first derivatives within a two-dimensional 
region R bounded by perimeter P and with first derivatives that are bounded at all 
point within R 

f/(r')V 2 g(r')rf 2 r' = 

(6) 

20 where n is the unit vector normal to P. If the function g(r f ) is replaced by the 
function G(r,r') satisfying 

V 2 G(r,r') = <y(r-r'), (7) 
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the LHS of Equation (6) becomes a formula for^r) in terms of its Laplacian inside R 
and boundary conditions on P: 

/(r) = £G(r,OVY(r^ (8) 

The function G(r,r') is called a Green's function and is defined by its Laplacian 
5 according to Equation (7). This means we are free to choose the boundary conditions 
for the Green's function. In general this choice is used to simplify Equation (8) by 
eliminating one of the two perimeter integrals. 

If the problem involves Neumann boundary conditions, 

n.V/(r) = A(r) for reP, 

1 0 where h(r) is known, the middle term of Equation (8) can be eliminated by specifying 

VG(r,r')n = 0 for reP. (9) 

The Green's function cannot simultaneously satisfy Equation (9) and Equation (7), as 
can be seen by applying Gauss' theorem: 

I V 2 G(r, r') d 2 r = <£vG(r, r') • fi dr . 

15 The RHS of this equation is identically zero, from Equation (9). The LHS must 
therefore also be zero. This is clearly not the case if Equation (7) is satisfied exactly; 
so instead we must have 

V 2 G(r,r') = <5(r-r')-i 
A 

where A is the area of the region R. This has the effect of subtracting the mean value 
20 of fir) on the LHS of Equation (8), which is acceptable since J[r) can only be 
determined to within an additive constant under Neumann boundary conditions. 

Referring to Equation (5) we see that the axial intensity derivative contains 
information about the normal wavefront slope around the aperture edge, rather than 
the function itself. The problem therefore involves Neumann boundary conditions in 
25 a natural way, so we adopt a solution with a Green's function satisfying Equation (9). 
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Inserting this into Equation (8) and replacing /(r) by <|>(r) gives the final solution for 
the wavefront phase as 

^(r)=| i G(r,r')VV(r')^V'-4 /) C?(r,r')V ( *(r')-nrfr' (10) 

Let us suppose that we can make a measurement of 

5 *(r>-|f/ z (r% (11) 

I Q dz 

If we multiply S(r') by G(r,r f ) and integrate over all space, we obtain (using Equation 
(5)): 

f £(r')G(r, r') d V = f [W A V V - S c n . V$>) G(r , r') rf V 

The aperture function ^ restricts the integral over VV to the region .4, and the 
10 delta-function 5 C converts the integral over n-V0 into a line integral around C. 
Therefore we have: 

jS(r')G(r, r') d V = £ G(r, r') V V(r') d V - <£ G(r, r ')V0(r') • n dr' (12) 

It can be seen that if we equate the regions A and R and the curves C and P, the RHS 
of Equation (10) is the same as the RHS of Equation (12). Thus 
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<*(r) = |5(r')G(r,r')rf 2 r' (13) 

gives the complete solution for the wavefront from a single area integral. 

As noted previously, many applications do not require the full wavefront #r), but one 
or more coefficients describing the amplitude of specific components, e.g. the amount 
of a particular Zernike mode or a set corresponding to the available modes of a 
2 0 wavefront modulator. If the wavefront is expanded as a linear combination of 
orthonormal functions «j(r), a particular coefficient is obtained by integration: 

a^\m<r)d 2 r (14) 
Substituting for $r) from Equation (13) gives 
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a x = \\S{r^G{x > r t )u i {r)d l r 1 d 2 r . (15) 

Defining: 

0,(0* |G(r,r'K(r)iV (16) 

Equation (15) becomes: 
5 a,= J^rOG^rOrfV (17) 

Therefore, the coefficient is obtained by the integral of the signal S(r) multiplied by a 
'modal projector function' G\{t\ 

Frequently the wavefront will be expressed as a vector, each value corresponding to a 
sample of the wavefront at a particular point, or the mean value of the wavefront over 
10 a small area. For sampling at a set of points {R,}, the corresponding expansion 
functions are: 

^.(rWtr-R,), 

and for mean values over a set of pixels centred on a set of points {Rj}, 

^(rJ^^r-R,)®^), 

15 where p(r) is the support function for one pixel centred at the origin. Once the set of 
expansion functions wj(r) is defined, the coiresponding set of modal projector 
functions can be calculated from Equation (16), and the a\ obtained by Equation (17). 

The function S(r') is approximated by the difference between two intensity 
distributions on either side of the aperture plane, separated by a small distance 8z: 

20 dl(r,z) ^ I(r,z + Sz/2)-I(r 9 z-Sz/2) 

dz Sz 

If the intensity distributions are to be recorded on a pixellated detector, the data will 
be in the form of a vector of intensity values, Sj, each being the mean value of Sir 1 ) 
over a pixel centred at r' = pj. The complete function is therefore approximated by 
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J 

where #(r ! ) is the support function for one pixel centred at the origin. Substituting this 
expression for S(r') into Equation (17) gives 



^/ = j[Z^( 5 ( r, -Py)®?( r ')f 



G,(rV r'. (19) 



5 Reversing the order of summation and integration and defining: 

G IJ = \(8{T'-p J )® q (T , ))G i {T')d 2 r' 

Equation (19) becomes 

0r-E*A/' (20) 
J 

which is a simple matrix-vector multiplication. The simplicity of this solution, 
10 combined with the fact that the modal projector matrix Gy can be pre-calculated, 
provides for the possibility of very fast modal wavefront sensing. 

The fact that the longitudinal intensity derivative is approximated by a difference in 
measured intensity distributions located a non-zero distance from the pupil plane 
means that there will be, in the measurement planes, non-zero intensity in the region 

15 exterior to the projected aperture. Thus, the modal projector function Gj(r) must be 
suitably defined in the region outside the aperture as well as within it. It has been 
found that a suitable scheme for circular apertures is to continue the boundary 
condition of Equation (9) to infinity. That is, the value of the Green's function at any 
point outside the aperture is taken to be the same as the value at the edge of the 

20 aperture at the same angle. Various other schemes were tested and this one was found 
to give the best performance in terms of minimising the error when sensing the low- 
order Zernike modes. 

Thus it will be seen that the re-formulation of the phase-diversity algorithm provides 
fast wavefront reconstructions and requires only a simple matrix multiply for data 
25 inversion. As implemented in experimental tests, the data collection is achieved in a 
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single image from a single, pixellated focal plane by use of a distorted diffraction 
grating. 

The algorithm is capable of working with both point sources and with extended, low- 
contrast scenes and can provide wavefront reconstruction data in any required set of 
5 basis functions. In particular, the formulation for point by point or for Zernike- 
polynomial decompositions is equally easily achieved. 

The algorithm can be used on wavefronts with severe aberrations, with extended 
sources, with partially obscured wavefronts, with strongly-scintillated wavefronts and 
with some cases where the wavefront is discontinuous (e.g. in multiply-connected 

10 pupils). It has already been tested in computer simulation and experimentally and has 
been found to be robust to many sources of experimental error. In most cases it is 
found experimentally that a departure from the theoretical restrictions, implied 
expressly or implicitly by the theory, leads to a wavefront reconstruction that is low- 
pass filtered. In this case the reconstruction of the lowest order modes, which are 

15 generally the most important sources of image degradation in the terrestrial imaging 
applications for which the algorithm was formulated, are found to be very accurate. 

When the modes are reconstructed serially on a 450 MHz PC, the data reduction takes 
50 n.sec per wavefront mode reconstructed. The best experimental results to date 
show a relative accuracy of better than 2 nm in the defocus terms when measured on a 
2 0 bright, monochromatic point source at 633 nm. 

As thus far described, it is found that when the wavefront sensor is used with a 
compact high brightness source on a dark background, such as a laser beacon imaged 
through a narrow-band filter, the accuracy of the results is good. 

When such a beacon cannot be used, a natural beacon, such as a sun glint or a bright 
25 object on a dark background is sought. However, in such circumstances the source 
power is unlikely to be concentrated in a narrow waveband, making it difficult or 
impossible to eliminate background radiation from the rest of the viewed scene, and 
this tends to lead to a loss in accuracy. This loss arises because the basic wavefront 
curvature signal is given by the difference in the two intensity profiles measured or 
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recorded either side of the system pupil, divided by their integration sum for 
normalisation: 



/,(r)-/ a (r) 
J(/,(r) + / 2 (r)) 



While a compact source (i.e. one completely contained within the field of view of the 
5 system) will provide intensity profiles Ii and h which lead to a correct curvature 
signal and hence correct wavefront mode coefficients, generally uniform background 
illumination will contribute identically to Ii and I2, thus contributing nothing to the 
difference signal (numerator) but increasing the normalisation factor (divisor). The 
result is that when the system so far described views a scene with a compact source 
10 on a uniform background it will provide curvature signals which have the correct 
spatial structure but which are inappropriately scaled, resulting in incorrect wavefront 
mode coefficients. 

By splitting the intensity profiles into a contribution from the compact source and the 
background: 

15 /,(r) = / w (r)+/ fl (r) , / 2 (r)=/ s>2 (r)+7 a (r) 

the curvature signal takes the form: 

S - BW "2j/ s (r) + 2j/,(r) 

where S and B indicate source and background contributions to the intensity profiles. 
By rescaling according to : 

f/,(r) + f/„(r) 

the correct curvature signal can be obtained. 

It follows that correct measurements can be made provided Is and I B can be 
measured. This is extremely hard to do in the pupil plane where measurements of Ii 
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and I2 are made. Nevertheless, since in a lossless system by the conservation of 
energy the integrated intensity in any plane perpendicular to the optic axis is invariant, 
the quantities Is and Ib can be measured in the image plane where source and 
background can be more easily separated. 

5 In a numerically simulated example, the scene consisted of a uniformly bright 
background with an off- axis source of twice the background intensity. A defocus of 
0,5 was assigned to all points within the scene and the wavefront sensor intensity 
profiles calculated as well as the image, which is shown in Figure 6 as an out-of focus 
spot 29 on a darker background 30. The defocus measurement without background 
1 0 correction was 0.07 waves, significantly in error. By applying background corrections 
as outlined above, the results were: 

Background level = 1 

Background integrated over unit circle: B = n 
Total integrated image: T = 3 .65 1 
15 Source S = T - B = 0.509 

Rescaling factor a = T/S = 3.651/0.509 - 7.17 

Corrected defocus measurement = 0.07 x a = 0.7 x 7.17 = 0.5, which is as expected. 

The reader is directed to our copending UK Patent Application No. GB 0301923.9 
(ref: P2 1386GB) which describes and claims alternative apparatus for providing the 
20 two focussed near- pupil- plane images in a common focal plane. The reader is also 
directed to our copending UK Application No. GB 0206242.1 (Ref: P2 1567GB) of 
even date with the present application which deals with an alternative method for 
implementing the matrix function by the use of a radiation (e.g. optical) grey-scale 
mask. 
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CLAIMS 

1. Measuring apparatus for determining data relating to the local shape (or 
distribution of local phase) of a radiation wavefront arriving at a pupil plane, 
wherein said shape is defined by a set of predetermined orthonormal 

5 functions, each function being provided with a weighting coefficient for 

determining the shape, said data comprising at least one said weighting 
coefficient, the apparatus comprising a said input pupil, rate means responsive 
to said radiation for determining a pixelwise distribution indicative of rate of 
radiation intensity change as the radiation traverses the input pupil, and 

10 converting means for converting said intensity distribution to said data, 

wherein said converting means comprises a store holding one or more 
matrices of predetermined values, each said matrix corresponding to one said 
orthonormal function, and the size of each said matrix corresponding to the 
number of pixels in said pixelwise distribution, and calculating means for 

15 multiplying said pixelwise distribution by a said matrix and adding the results 

to provide said weighting coefficient for its said orthonormal function. 

2. Apparatus according to claim 1 wherein said rate means comprises optical 
means for providing focussed first and second images of respective first and 
second secondary planes lying adjacent, and respectively before and after, the 

20 pupil plane, and deriving said pixelwise distribution as the pixelwise 

distribution of intensity difference between the first and second images. 

3. Apparatus according to claim 2 wherein one said secondary plane is a virtual 
plane. 

4. Apparatus according to any preceding claim wherein the optical means is 
2 5 arranged to provide said first and second images in a single plane. 

5. Apparatus according to any preceding claim wherein the optical means 
comprises a distorted diffraction grating or a computer generated hologram. 

6. Apparatus according to any preceding claim wherein the set of orthonormal 
functions are Zernike polynomials. 
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7. Apparatus according to any preceding claim adapted for use with optical 
radiation. 

8. Apparatus according to any preceding claim wherein said store holds more 
than one said matrix. 

5 9. Apparatus according to claim 8 wherein said store holds a number of matrices 
less than the number of pixels in said pixelwise distribution. 

10. Apparatus according to any preceding claim, and according to claim 4, and 
further including illuminating means for illuminating a point in the field of 
view, the apparatus being arranged so that light from said point in said 

1 0 secondary planes is spatially separated in said single plane. 

11. Apparatus according to claim 10 wherein the illuminating means is arranged 
to illuminate a plurality of said points in the field of view, the light from each 
point in each of said secondary planes being spatially separated in said single 
plane. 

15 12. Apparatus according to any preceding plane and including means for 
measuring the degree of focus, or range, of an object in the field of view as a 
function of at least one said weighting coefficient. 

13. Apparatus according to any preceding claim and further comprising a spatial 
phase modulator for receiving radiation from said input pupil, and modulator 

20 control means responsive to the said weighting coefficient(s) for controlling 

the modulator. 

14. Apparatus according to claim 13 wherein the spatial modulator is located 
between the input pupil and said rate means. 

15. Apparatus according to claim 14 wherein said modulator control means is 
2 5 arranged to minimise said weighting coefficients. 

16. Apparatus according to any one of claims 13 to 15 and further comprising 
imaging means for receiving radiation from said modulator. 
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17. Apparatus according to claim 13 and further comprising imaging means for 
receiving radiation from said modulator, wherein said rate means is arranged 
to receive radiation which has not been transmitted by said modulator. 

18. Apparatus according to any one of claims 13 to 17 wherein said modulator 
includes a tip/tilt mirror. 

19. Apparatus according to any one of claims 13 to 18 wherein said modulator 
includes a deformable mirror or a deflectable mirror array, or a liquid crystal 
device. 

20. Apparatus according to any preceding claim for use with an image comprising 
at least one relatively bright small area or spot on a background, the apparatus 
comprising means for measuring the brightness of the bright area or spot and 
the brightness of the background, and for combining the two brightness 
measurements to derive a scaling factor for application to a measured 
weighting coefficient. 

21. Apparatus according to claim 20 wherein the brightness measurements are 
effected in an image plane. 

22. A measurement method for determining data relating to the local shape (or 
distribution of local phase) of a radiation wavefront arriving at a pupil plane, 
wherein said shape is defined by a set of predetermined orthonormal 
functions, each function being provided with a weighting coefficient for 
deteraiining the shape, said data comprising at least one said weighting 
coefficient, the method comprising deteraiining a pixelwise distribution 
indicative of rate of radiation intensity change as the radiation traverses the 
input pupil, and converting said pixelwise distribution to said data, wherein 
said converting step comprises providing one or more matrices of 
predetermined values, each said matrix corresponding to one said orthonormal 
function, and the size of each said matrix corresponding to the number of 
pixels in each of said first and second images, for each said matrix multiplying 
said pixelwise distribution by said matrix and adding the result to provide said 
weighting coefficient for its said orthonormal function. 
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23. A method according to any one of claim 22 wherein the radiation is optical 
radiation. 

24. A method according to claim 22 or claim 23 and including the step of 
illuminating or irradiating a point in the field of view, and determining said 

5 data from radiation received from said point. 

25. A method according to claim 24 wherein a plurality of separated points in the 
field of view are illuminated or irradiated, and said data is determined from 
radiation received from each said point. 

26. A method according to any one of claims 22 to 25 and comprising the further 
10 step of modulating radiation from said input pupil with a spatial phase 

modulator controlled in response to the said weighting coefficient(s). 

27. A method according to claim 26 wherein the step of determining said 
pixelwise distribution is performed on radiation transmitted by said spatial 
modulator. 
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